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Abstract 

The main goal of the paper is to introduce methods which compute Bezier curves 
faster than Casteljau's method does. These methods are based on the spectral fac- 
torization of a n X n Bernstein matrix, B^{s) = PnGn{s)Pn^ , where P„ is the nxn 
lower triangular Pascal matrix. So we first calculate the exact optimum positive 
value t in order to transform P^ in a scaled Toeplitz matrix, which is a problem 
that was partially solved by X. Wang and J. Zhou (2006). Then fast Pascal matrix- 
vector multiplications and strategies of polynomial evaluation are put together to 
compute Bezier curves. Nevertheless, when n increases, more precise Pascal matrix- 
vector multiplications allied to affine transformations of the vectors of coordinates 
of the control points of the curve are then necessary to stabilize all the computation. 
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1 Introduction 



Bezier has his name on the curve B of degree n — 1 defined from n given points 
Zq = (xo,yo), Zi = Zn-i = in as follows: 

n-l 

B{s) = {x{s),y{s)) = ^ ZA,n-i(s), s e [0, 1], 

1=0 
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where 6j_„_i(s) = — s)" ^ for each i G {0, — which is called 

a Bernstein polynomial P]. That is, for each s G [0, 1] we have 

n— 1 71—1 

^(s) = bi^n-i{s)xi, y{s) = &i,„_i(s)?/i. 

Notice that B{s) = Br{l — s), where Br{s) is the Bezier curve defined from 
the control points in the reverse order: 
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Let s be a real number. The n x n lower triangular matrix B^{s) such that 
{B^)^j (s) = j_i(s), for each n > i > j > 1, is called a Bernstein matrix 

Paul de Casteljau developed a very stable algorithm to evaluate Bezier curves. 
In this so called Casteljau's algorithm the number of arithmetic operations 
grows quadratically with n: it requires n{n — l)/2 additions and n{n — 1) 
multiplications to calculate a point (which is not an endpoint) on a Bezier 
curve of degree n — 1 [3]. A natural question is if there could be a less expensive 
algorithm to compute a Bezier curve, which is answered e.g. in [5] for n < 10. 
Here, we extend this answer to n < 64 by a different approach, which arises 
from the expression of a n x n Bernstein matrix in terms of the lower triangular 
Pascal matrix P„ (see pQ) as follows: 

x{s) = elPnG{-s)PnG{-l)x, y{s) = e^P„G'(-s)P„G'(-l)y. 

Here, x = (xq Xi ■ ■ ■ y = {yoUi ■ ■ ■ 2/n-i)^; e„ denotes the nth canonical 

vector of M"; G{s) is the diagonal matrix such that e^G{s)ek = s^~^ for all 
k G {1, .■■,n}; and P„ is defined by 



{Pn 



, for ^ ^ j; 
, otherwise. 



From now on, the methods introduced here which are originated from the 
above expression will be called Pascal matrix methods. One central problem 
in these methods is how we can accurately compute a matrix-vector multi- 
plication with the n X n lower triangular Pascal matrix. In §2, we present an 
algorithm that utilizes only n(n-l)/2 additions that computes a matrix-vector 
multiplication with P„ in a very precise way, which is based on the property 
of the lower triangular Pascal matrix be a product of bidiagonal matrices. But 
this time complexity can be reduced to C(nlog n) by using the fact that the 
lower triangular Pascal matrix is similar to lower triangular Toeplitz matri- 
ces via diagonal matrices [B]. This transformation depends on a positive real 
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parameter t, which is arbitrary. In [6], it is suggested that {n — l)/e could 
be taken as a good approximate value for the optimum parameter. Here, the 
optimum value when it exists is calculated. It is left as an open problem if the 
set of n for which it does not exist the optimum value is finite. At the end 
of the section, there are results of tests to compare the accuracy of matrix- 
vector multiplications carried out by the fast multiplication algorithm, with 
both approximate and exact optimum values, and by the algorithm based on 
the bidiagonal factorization. 

The other central problem in Pascal matrix methods is, once calculated z = 
PnG{—l)x (or z = PnG{—l)y), how to efficiently evaluate e^PnG{—s)z. In 
§3, some results obtained with the use of the 0{n\og n) Pascal matrix- vector 
multiplication algorithm coupled with a Horner-type scheme in the computa- 
tion of Bezier curves of degree n — 1 are then presented. When n increases, the 
fast matrix-vector multiplication becomes unstable, as well as the evaluation 
of B{s) for s close to 1. Then, the matrix-vector multiplication must be done 
in a more precise way, and a stabilizing procedure to -B(s)-evaluation should 
be attempted, e.g. by dividing the evaluation process in two: first, compute 
B{s) for s G [0, 1/2); then, compute Br{s) for s G [0, 1/2], where B^^s) is the 
reverse Bezier curve, that is, Bj.{s) is the Bezier curve defined from the points 
Pn-i, -Pi and Pq, in this order. For n > 32, even the 2-steps -B(s)-evaluation 
becomes unstable, yielding incorrect values for s close to 1/2. It is when we 
introduce an affine transform of the vectors of coordinates in order to improve 
the evaluation. 



2 Pascal matrix-vector multiplication 



In this section some algorithms of Pascal-type matrix-vector multiplication 
are discussed from their time complexity. These algorithms are founded on 
algebraic properties of these matrices which are also presented in the following. 



2.1 (9(n^) arithmetic operations methods 



Algorithms of matrix- vector multiplications with P„, PnGnit) and B^{t), re- 
spectively, are here presented in the form of MATLAB functions, all of them 
demanding 0{n'^) arithmetic operations. 
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We begin observing that 
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that is 
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E2-E\ — E^.E2-Ei. 



Since for 1 < j < i < n 



we conclude that 



1 

OPn-l 



i - 1 



1 

An-l) 



Eji-i 



Therefore, it has just proved by induction the following statement: 

Proposition 1 Let P„ be the n x n lower triangular Pascal matrix. Then 
Pn = En-i El where, for all k G {l,...,n-l}, E^ = I+ek+iel+...+ene^_^. 



The first algorithm is described as a MATLAB function in the following: 
function x = pascal_product (x) 

y.PASCAL.PRODUCT Multiply a vector x by the lower triangular 
"/oPascal matrix 
n= length (x) ; 
X = x( : ) ; 
for k = 2:n 

for s = n: -1 :k 

x(s) = x(s) + x(s-l) ; 

end 

end 
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Notice that from the above factorization we also conclude that 

Ik \fPkO 

Pn^k \ Z W 



where 



Z W 



" (9 - ( 



\ .. 

A similar factorization has PnGnit): 



k\ (k 



0/ VI 



A 





1}J 



( 



1 

1 
1 G)^ G)^^ 









^ ... 0^ 
1 ... 



1^1 (v)t t.y ... (::;)t"-v 

that is 



0... 10 
0... It 



1 ... 0^ 
1 t ... 
1 t ... 



0... It 



Proposition 2 Let P„ he the nxn lower triangular Pascal matrix and G{t) = 
diag{l,t,t^,...X~^)- Then PnGn{t) = En-i{t)...Ei{t) where, for 1 < k < 
n — 1, 



Ek{t) = eief + ... + Ckcl + ek+i[ek + tck+i]'^ + ... + e„[e„_i + te„]'^. 



The second algorithm, which utilizes Proposition [21 is displayed just below: 
function x = pascal_g_product(x,t) 

7oPASCAL_G_PRODUCT Multiply a vector x by PG(t) where 

°/o P is the lower triangular Pascal matrix and 

7. G(t) = diag(l,t,t~2, . . . ,t~{n-l}) 

n= length (x) ; 

X = x( : ) ; 

for k = 2:n 

for s = n:-l :k 

x(s) = x(s-l) + t*x(s) ; 

end 

end 
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By using the notation of [4J, a Bernstein matrix can be described as 

B:{t)=P„\l-t]Gn{t), 

where {Pn[t])ij = {Pnliji^^^i 1 < j < * < Therefore, it is not difficult to 
conclude that a Bernstein matrix has the following bidiagonal factorization: 



1 

1 -t 

















n-3+2 (n-l\^n-l 



\ 



10... 
1... 







1 



... l-t t 



1 ... 
l-t t ... 
l-t t ... 



V 



...l-tt 



that is 

Proposition 3 Let B'^{t) be a n x n Bernstein matrix. Then 

B^{t) = El_^{t)...El{t) where, for I < k < n - I, 
El{t) = cicl + ... + ekel + ek+i[{l-t)ek + tek+if + ... + en[{l-t)en-i + ter. 



Proposition [3] yields the following algorithm, that is essentially the Casteljau's. 
function x = bernstein_product (x,t) 

7oBERNSTEIN_PRODUCT Multiply a vector x by a Bernstein matrix 
n= length (x) ; 
X = x( : ) ; 
tl = l-t; 
for k = 2:n 

for s = n:-l :k 

x(s) = tl*x(s-l) + t*x(s); 

end 

end 
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Similarly we can conclude that 



\f Blit) 

^ ] \ z w 



) 



where ( Z I is the matrix 



{i-tf (9(1 



tf-H 
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2.2 Balancing a lower triangular Pascal matrix 



Let Zq = {xo,yo), Zi = {xi,yi), Z„„i = (x„_i,y„_i) be n points of M^. 
Let X = (xqXi ... Xn-iV and y = {y^y^ . . . y^-i)^ . For s G [0, 1], let B^{s) 
be the n x n Bernstein matrix. A spectral decomposition of B^{s) is B^{s) = 
PnGn{s)P~^, where P„ is the nxn lower triangular Pascal matrix and = 
diag{l, s, s"-^) ([1]). Since P'^ = G„(-1)P„G„(-1) ([1]), we conclude that 
the coordinates of the Bezier curve B{s), s G [0, 1], defined from Zq, Zn-i 
(in this order) is given by 



Xis) = elPnGn{~s)PnGn{-l)x, y{s) = elP^G n{- s) PnGn{-l)y . (1) 



The computation of a matrix-vector multiplication with the Pascal matrix 
Pn can be carried out in 0{n\ogn) operations |6]. Hence, for each s G [0, 1], 
the computation of x{s) and y{s) also requires C(nlog?T,) operations. How- 
ever, arithmetic operations with Pascal matrices are very unstable because of 
the various magnitudes their entries have. In this section, a preconditioning 
technique is proposed to deal with this instability. This technique, which is 
introduced in [6], is based on the factorization 



Pn = Dr,{t)Tn{t) Dnit) 



(2) 
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where Dn{t) = diag{0\ 



1! (ri-1)! 



fn— 1 



), and 



Tn{t) 



1 
t 

V. 




1 








^n— 1 jn — 2 ^n — 3 

V (n-1)! (n-2)! (n-3)! 



...o\ 

... 
... 



Note that for < t < 1 the entries of the Toephtz matrix T„(t) vary from 1 
to that is to say T„(t) is ill-conditioned. In [6] it was found that a good 



n 



(n-l)!' 

value for t is ti 

e 

almost all values of n: to 



Here we propose a more accurate value for t, for 
where 



'(n-l)! 



1)!J. 



In order not to cause great instability in matrix-vector multiplication with 
Pascal matrices, we will balance these matrices by using the decomposition 
given in ([2]). Since the entries of the matrix Tn{t) are of the form 

f{m) = — , 
ml 

where m = 0, 1, ... ,n — 1, we would like to find out a value t for which 
max / / min / is the least possible. If this optimum t exists, it will bring those 
entries to be the closest in magnitude to one another. The following lemma 
resumes basic facts about the monotony of /. 

Lemma 4 Let t G (0, oo). Consider the function / : N — ?■ [0, oo) defined by 

f{m) = — -. Then 
ml 

(i) f is a nondecreasing function for integers m such that m ^ \t\ ^ t; 
(a) f is a nonincreasing function for integers m such that t ^ \t~\ ^ m. 



PROOF. 

(i) f{m - 1) 

(ii) /(m-1) 
□ 



m—l 



m t 



m—l 



(m-1)! 



m (m 



1 ! 



m f' 



(m — l)! m(m— 1)! ml 



f{m), since m ^t. 
f{m), for m ^ t. 



8 



Therefore, if the optimum value of t exists it would belong to the interval 
[1, n — 1). Now, max // min / = , ii < m < t; if t < m < n — 1, then 

( L'- J ) • 

t\^^(n - 1)! 

max f/ min f — ^ , ^ -, , , . That is, for the optimum value 

is achieved. The following lemma is about the monotony of /i and f2- The 
proof is an exercise of Calculus and will be omitted here. 

Lemma 5 Let n he an integer greater than 1. Let fi and /2 be two functions 

defined at the subset (0,n — 1) of real numbers such that fi{t) — JJtJy. 

t^'Xn-l)\ , 
f2{t) = -—qj^, respectively. Then 

(i) fi is an increasing function; 
(a) /2 is a nonincreasing function. 

Moreover, fi is continuous and /2 is continuous but at integer values. 

Remark 6 Therefore, if exists f e (0,n — 1) such that fi{i) = f2ii), this is 
the optimum value. 

Now, we can calculate the optimum value when it exists. 

Lemma 7 // exists t e (A;, A; + 1) such that = f2it), then 



V 

PROOF. Let t be such that fi{t) = /2(f). Since i e (A;, A; + 1), it follows that 



Ik k + l 1 ~ „-2/(n-l) 



k\ (A; + l)!t"-i (n-1)! t«-2 V A; + 1 

□ 

Lemma 8 // there is a i e (A;,A; + 1) such that fi{i) = f2{i), then k — 

(n-lV 

PROOF. Prom the former lemma, if A; < t < A; + 1 then = -. 

k + l 
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Hence, 

A;-2<l!pl):<(A; + l)"-2. 
/b ~r -L 

But this yields 

{n - 1)! <{k + 1)"-^ and (n - 1)! > {k + l)A;"-2 > 
Therefore, 

k = max{m e Z+| m'^"^ < (n - 1)!} = max{m e Z+| m < ""-{/(n- 1)!} 
That is. A; = L ""^(n- 1)!J . □ 



Remark 9 Note that, for n = 3, = 1, /2(1) = 2, /2(1+) = 1. So, there 

is not the optimum value because /2(1) > /i(l) > h{^~^) (^nd a number greater 
than 1 but very close to it can be taken as a good value to balance the 3x3 

Pascal matrix. There are several positive integers n, n > 3, such that for some 
integer k, 1 < k < n — 2, f2{k) > fi{k) > f2{k~^). And this happens if and 
only if 

k''(n-l)\ k^ ^ k''+\n-l)\ 
k\ A;"-i ^ ¥ - (A; + l)!A;"-i' 

that is, 

< (n - 1)! and k""-^ + k""'^ > {n - 1)!. 



Proposition 10 // there is some integer number n, n > 3, such that ^ < 
(n — 1)! and k"'~^ + > (n — 1)! for some integer k, 1 < k < n — 2, then 

n — 1 , n — 1 
< k < — - — . 



PROOF. If A; > (n - l)/2, then 



Since, for n > 1, n** > (n!)^, we have 



„_i^,n-L (n-l)-! 



(n - ^ ((r? - 1)!)^ 



2n— 1 — 2"'"''" 

Now, for n > 3, 2"-^ < (n - 1)!. Hence k^-^ > {n - 1)!. 
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On the other hand, k must be greater than (n — l)/e, because 



77,-1 



< V'2vr(n-l)(^)"-i<(n-l)!, 



by the Stirhng formula. □ 



160 



140 - 



120 



100 - 




5.5 5.6 5.7 5.8 5.9 6 6.1 6.2 6.3 6.4 6.5 



Fig. 1. n = 15: /i and /2 don't intersect 

Remark 11 As far as we know, the conjecture about the set of integers n 
such that k"'~^ < (n — 1)! < k"'~^ + fc""^ for some integer k, 1 < k < n — 2, 
be finite has not been proved yet. Each integer n less than 10000 belonging to 
this set is displayed together with its corresponding integer k, 1 < k < n — 2, 
in TableUl 

For each value of n in the table, we will take t = k as the optimum value. For 
the other values of n, 4 < n < 10000, the optimum value is given by 



(n-1)! 
k+1 



where k = [" \J {n — 1)!J . 



In Table [2] we can compare the accuracy of the computation of Pz for some 
values of n, where z is the vector defined by Zk = (—1)*^ for k = l,...,n: 
wi = Dn{l).Tn{l).Dn{l)-'z; W2 = D^{t) .Tn{t) .D^{t)-^ z , with t = {n- l)/e; 
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n 


k 


iO 


D 


oy 


iO 


7/1 


Zo 


OZ ( 


1 Qi^ 

iyo 


"^l 71 


i iuo 




91 71^ 

Z 1 / (J 


7036 


2590 


7534 


3194 


7537 


3401 



Table 1 

^n-i < _ 1)1 and /c"-! + > (n - 1)! 

ti;3 = Dn{t).Tn{t).Dn{t)-^z, with t = ""^^^^ and = [ "-^(n-l)!J, or 

t = A; if (n, /c) is one of the pairs in Table [11^4= pascal jproduct{z). For 
the computation of Wi, W2 and W3 a fast Toephtz matrix- vector muhiplication 
was used, which demanded only 0{n\ogn) operations. Remark that in exact 
arithmetic PnZ = ei, the first canonical vector. 



n 


\\PnZ -Wi\\2 


\\PnZ - W2II2 


\\PnZ -W3\\2 


\\PnZ - W4\\2 


4 


4.9262e-16 


4.4790e-16 


9.6908e-17 





8 


2.4685e-13 


9.5720e-15 


6.2521e-15 





15 


1.7782e-06 


5.4617e-13 


9.3638e-13 





16 


1.9051e-05 


9.5411e-13 


8.4432e-13 





32 


1.0318e+16 


2.7147e-07 


9.5906e-08 





39 


2.8076e+27 


1.6654e-05 


4.4418e-06 





64 


2.0661e+71 


1.7187e+04 


7.2813e+03 






Table 2 

Errors in Pascal matrix-vector multiplication 



3 On computing a Bezier curve 

In this section we are going to compute a Bezier curve. First, notice that if the 
control points are translated by a vector v = {p,q), the Bezier curve is also 
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translated by v. Another feature of a Bezier curve is that an uniform scahng 
of the control points yields an uniform scaling of the curve. Hence, without 
loss of generality we assume that the coordinates of the control points are all 
positive and less than or equal to 1 and, therefore, the oo-norms of P^x and 
PnU are less than or equal to 2". 

3.1 B{s)- evaluation 

Pascal matrix methods compute a Bezier curve B{s) of degree n — 1 via the de- 
composition -B^(s) = PnGn{—s)PnGn{—^)'- first z = P„G'„(— l)x = PnX- and 
w = PnGn{—l)y = PnU- are computed; then the polynomials e^PnGn{—s)z 
and e^PnGn{—s)w are evaluated. When n is small, e.g. n = 32, z and w 
have 2-norms around 10^ for ||x||oo = 1 or ||y||oo = 1? and both polynomials 
could be efficiently evaluated for each s G [0, 1]. The function g_pascaLproduct 
does the evaluation, with n{n — l)/2 additions and n{n — l)/2 multiplications 
for each s. A less expensive alternative is to use a Horner-like scheme that 
evaluates the polynomial concomitantly with the binomial coefficients. 



The test control points were defined by the MATLAB function rand{n,2), 
which returns an n x 2 matrix containing pseudo-random values drawn from 
a uniform distribution on the unit interval. We have used the Pascal matrix- 
vector multiplication done from the similar Toeplitz matrix T{t), with t found 
by our procedure (according to our calculations in the last section), plus the 
B{s) evaluation given by the Horner-like scheme cited above. The results were 
compared with the ones obtained by Casteljau's method. In table [3] the cxd- 
norms of the differences of the results obtained by the Pascal matrix method 
(Bp) and Casteljau's (Be) are listed close by the average time of computation 
of all B(s) evaluation {As = 1/128). We have taken the smallest time among 



n 


\\Bp — -Be oo 


timep 


timec 


4 


7.7716e-16 


0.001810 


0.001273 


8 


2.8547e-14 


0.001954 


0.002164 


15 


9.3585e-ll 


0.002095 


0.005024 


16 


1.9592e-10 


0.002238 


0.005475 


24 


1.2341e-06 


0.002518 


0.009927 


32 


0.0190 


0.003169 


0.015896 



Table 3 

Pascal matrix method x Casteljau's 
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10 elapsed times obtained from consecutive executions of the procedure as the 
average time of computation, all of them computed by the MATLAB's built-in 
tic/toe functions. 

For n = 32, the i?(s)-evaluation becomes unstable when s approaches to 
1. To locally and globally improve the evaluation, we have made a simple 
procedure which has yielded more precise results, that is to divide the process 
of evaluation in two independent steps: 

(a) evaluate e^P„G„(— s)z and e^P„G'„(— s)w for < s < 1/2; 

(b) evaluate e^PnGn{—s)zr and e^PnGn{—s)wr for 1/2 > s > 0, which is 
equivalent to evaluate e^PnGn{—s)z and e^P„G„(— for 1/2 < s < 1. 

The results indicate that the procedure has worked well as far as n = 41 and 
some of those can be seen in tabl^H 



n 


\\Bp — BcWoa 


timep 


timec 


32 


2.3113e-07 


0.005602 


0.015901 


36 


2.6961e-05 


0.006106 


0.019464 


39 


1.3152e-04 


0.005873 


0.022314 


41 


4.8668e-04 


0.006229 


0.024372 


42 


0.0022 


0.006387 


0.025476 


48 


0.1112 


0.006706 


0.032366 



Table 4 

Pascal matrix method with reverse evaluation x Casteljau's 

Two remarkable facts arise in table HI first, the time of computation for n = 39 
is smaller than that for n = 41; second, the sudden loss of precision from 
n = 41 to n = 42. One of the explanations for the first fact is because there is 
no calculation to find the optimum value to balance P39: 15 is taken to be the 
optimum value. The second fact surely has to do with the limitations of the 
flowing point arithmetic of our machine, a 32-bits AMD Athlon XP 1700+ 
(1467 MHz). 

3.2 On conditioning the vectors of coordinates 

One way to overcome this lack of stability is to transform the vectors x and y 
of coordinates into a vector very near to e^^ = (1 1 ... 1)"^. Since 

1 m 

Vk = rr^k-i H —re 

m + 1 m + 1 
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is a stationary scheme that converge to the solution e of /x = e for any vq, the 
idea is to compute the Bezier curve T^^B) from control points Wq = Tm(Zo), 
Wn-i = Tm{Zn-i), whcrc Tm{v) = {v + m.e)/{m + 1), and then to obtain 
B by inverse transforming the points of T{B). Note that 

Tr^{B){s) = elP^G{-s)Pr.G{-l) (t„(x)//T„(i/) ) • 

Hence, 

elPnG{-s)PnG{-l)T^{x) = e^P„G(-s)P„G(-l)(x + m.e)/(m + 1) = 

1 T?7 

= —-elP^G{-s)Pr,G{-l)x + —-elPr,G{-s)Pr,G{-l)e = 
m + 1 m + 1 



m + 1 m + 1 

Thus, since an analogous result is obtained with y, x{s) = (m + l)Tm(_B)(s) — 
m. Observe that 

TmTm'{v) = {v+ [{m + l)(m' + 1) - l].e) /[(m + l)(m' + 1)] = T(jr,+i){m'+i)-i- 
Our strategy to have better i?(s)-evaluation is the following: 

(a) evaluate e^PnGn{—s)z and e^P„Gn{,—s)w for < s < 1/3; 

(b) evaluate e^P„G„(-s)P„G'„(-l)T^(x) and e^P„G'„(-s)P„G„(-l)T„(?/) 
for 1/3 < s < 2/3, and inverse transform; 

(c) evaluate e^PnGn{—s)zr and e^PnGn{—s)wr for 1/3 > s > 0. 

For the experiments, the coordinates of the n control points were defined from 
the command A = rand{n, 2), followed by the normalization A = A/norm{A). 
Tm o Tm, where m + 1 = 32768 (2"'^^), were applied to A when n = 42, 48, 54; 
Tm' o Tm o Tm, whcre m' + 1 = 1024, when n = 59; Tm" o Tm' o Tm o Tm, where 
m" + 1 = 4, when n = 64. The results are in tabl^ 



n 


\\Bp — -Be oo 


timep 


timec 


42 


8.3290e-07 


0.013525 


0.025363 


48 


1.7620e-06 


0.015566 


0.032190 


54 


2.3903e-04 


0.017786 


0.039937 


59 


9.9235e-04 


0.019782 


0.047320 


64 


0.0048 


0.022044 


0.055178 



Table 5 

Pascal matrix method with piecewise evaluation x Casteljau's 
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4 Conclusions 

We have presented results obtained from some methods to compute a Bezier 
curve of degree n — 1, for various values of n. They were created from a 
description of the curve that involves matrix-vector multiplications with the 
n X n lower triangular Pascal matrix P„, which are here called Pascal matrix 
methods. With this in mind we have introduced two algorithms: one, which 
only demands n{n — l)/2 additions, is very precise and it is based on the fact 
that Pn is a product of bidiagonal matrices with and 1; the other, which 
demands 0{n log n) algebraic operations, depends on a positive real value 
in order to minimize the magnitudes of the entries of P„ when considered as 
a scaled Toeplitz matrix. We have seen that there is a function that relates 
n to that optimum value except for the integers belonging to a certain set, 
which we have not yet known if it is finite or not. Once the matrix- vector 
multiplication done, a polynomial evaluation should have been carried out for 
various values s G [0, 1], which has become unstable as s approaches 1/2. From 
the set of experiments presented here we have seen that the combination of 
Pascal matrix- vector multiplication plus polynomial evaluation has converged 
to the results obtained from Casteljau's by adopting some strategies, which 
vary according to the magnitude of n. And even so, they are more effective 
concerning time of computation than Casteljau's, at least for n < 60. 
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